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ABSTRACT 

We have investigated the stability of a set of non-rotating anisotropic spherical 
models with a phase-space distribution function of the Osipkov-Merritt type. The 
velocity distribution in these models is isotropic near the center and becomes radially 
anisotropic at large radii. They are special members of the family studied by Dehnen 
and Tremaine et al. where the mass density has a power-law cusp p oc r"^ at small 
radii and decays as p oc at large radii. 

The radial-orbit instability of models with 7 = 0, 1/2, 1, 3/2, and 2, was studied 
using an N-body code written by one of us and based on the 'self-consistent field' method 
developed by Hernquist and Ostriker. These simulations have allowed us to delineate a 
boundary in the (7, ra)-plane that separates the stable from the unstable models. This 
boundary is given by IT^jTi = 2.31 ±0.27, for the ratio of the total radial to tangential 
kinetic energy. We also found that the stability criterion df /dQ < 0, recently raised by 
Hjorth, gives lower values compared with our numerical results. 

The stability to radial modes of some Osipkov-Merritt 7-models which fail to satisfy 
the Doremus-Feix's criterion, df/dE < 0, has been studied using the same N-body 
code, but retaining only the / = terms in the potential expansion. We have found no 
signs of radial instabilities for these models. 



Subject headings: galaxies: kinematics and dynamics — instabilities — methods: nu- 
merical 
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1. Introduction 

Spherical stellar systems in which the orbits are strongly eccentric, or radial, are unstable to 
forming a bar. This so-called radial-orbit instability (hereafter ROI), was first demonstrated for a 
sphere consisting entirely of radial orbits by Antonov (1973). It was actually observed by Henon 
(1973) in a series of N-body simulations for the generalized polytropes, defined by a distribution 
function which is the product of power laws in energy and angular momentum. Subsequent N- 
body simulations (Merritt & Aguilar 1985; Barnes, Goodman, & Hut 1986; Dejonghc &: Merritt 
1988) and numerical linear stability analysis (Saha 1991; Weinberg 1991) have revealed that a small 
velocity anisotropy can be enough to make a system unstable; they also have showed that such a 
system evolves quickly into a triaxial bar. 

The ROI can be explained by the transformation of loop orbits into box orbits when the 
spherical symmetry of the potential is broken (Merritt 1987). The mechanism is similar to the 
scenario described by Lynden-Bell (1979) for the formation of a bar in a disk. Initially, all the 
orbits are precessing ellipses (loops). A weak bar-like (/ = 2) perturbation generates a small 
torque that could trap those orbits with the lowest angular momenta; realign them around the bar, 
reinforcing the strength of the initial perturbation. This instability implies an upper limit to the 
degree of velocity anisotropy in spherical models. Dejonghe k, Merritt (1988) studied the stability 
of the Plummer (1911) model with two different anisotropic distribution function (hereafter DF) 
and concluded that the ROI is the most robust and is likely to be the most important instability 
in elliptical galaxies. 

For spherical stellar systems with isotropic DF, / = f{E), it is known that df/dE < acts 
as a sufficient condition for the stability under radial and nonradial perturbations (Antonov 1962; 
Sygnet et al. 1984). However, the situation for anisotropic spherical models is poorly understood 
and only the stability to radial modes can be analytically tested using the sufficient condition 
df/dE < (Doremus & Feix 1973). A sufficient stability criterion based in the global anisotropy 
parameter 2Tr/Tt, the ratio of kinetic energies corresponding to the radial and tangential direc- 
tion, has been suggested by Polyachenko & Shukhman (1981; see Fridman Sz Polyachenko 1984). 
However, N-body simulations have showed that is not a reliable criterion because its value seems 
to be model dependent (Merritt & Aguilar 1985; Dejonghe &: Merritt 1988). Moreover, Palmer & 
Papaloizou (1987) have analytically demonstrated that no stability criterion can be found using 
the anisotropy parameter 2Tr/Tt, because there is always a spectrum of unstable (growing) modes 
if the DF becomes singular for a vanishing angular momentum; no matter how weak appears the 
DF divergence. 

Recently, Hjorth (1994) has proposed a simple criterion for the stability of models with DF 
of the Osipkov-Merritt type (Osipkov 1979; Merritt 1985). These models are described by an 
anisotropic DF, f^E^L?) = f{Q), with Q = E — L^/2r^; where E is the energy per unit mass 
and L is the magnitude of the angular momentum per unit mass. Inside the anisotropy radius ra 
the velocity distribution is nearly isotropic, while outside ra the radial anisotropy increases. Thus, 
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the parameter controls the global anisotropy degree in the velocity distribution; the smaller its 
value, the larger the anisotropy of the system. 

In this paper we show the results of a numerical study of the stability properties for a set of 
Osipkov-Merritt models. These are special members of the family of spherical 7-models indepen- 
dently studied by Dehnen (1993, hereafter D93) and Tremaine et al. (1994, hereafter T94); with 
densities that behave as p oc r~'^ near the center and as p oc r~'^ for large radii. The stability 
was investigated using an N-body code based on the 'self-consistent field' method developed by 
Hernquist k Ostriker (1992, hereafter H092). 

In § 2 we briefly describe the 7-models and the N-body code used to study their evolution. In 
§ 3 we present our results for models with 7 = 0, 1/2, 1, 3/2, and 2. In particular, we are able to 
establish a stability boundary for the ROI in the (7, ra)-plane. A discussion and summary of our 
results appear in § 4. 

2. N-body Simulations 
2.1. Models and Initial Conditions 



In this section, we briefly summarize the properties of the one-parameter family of spherical 
models studied by D93 and T94 (more details can be found in the references just quoted). These 
7-models have a mass density given by 

3-7 aM 



p{r) 



Att r')'(r-ha)4-7' 



< 7 < 3, 



(1) 



where a is a typical scale length and M is the total mass of the system (T94 describe the same 
density by the parameter r/ = 3 — 7 and call them "r/- models" ) . The models of Jaffe (1983) and 
Hernquist (1990, hereafter II90) are recovered for 7 = 2 and 7 = 1, respectively. 

The gravitational potential associated with the density ([l|) is obtained from the Poisson equa- 
tion, and has the simple form 
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where G is the gravitational constant. In the following, we adopt a system of units in which M, a, 
and G are unity. 

Many dynamical properties of the 7-models with an isotropic DF, / = f{E), were given by 
D93 and T94. The former also gave analytical expressions for anisotropic DFs of the Osipkov- 
Merritt type for some 7-models. More recently, Carollo, de Zeeuw, & van der Marel (1995) have 
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been able to express the projected line-of-sight velocity profiles of models with 7 < 2 in a single 
quadrature. Models with 7 > 2 are less suited to describe spherical galaxies because they display 
some unrealistic properties, such as infinite central potential and surface density profiles that differ 
significantly from the i?^/^ law (see D93). 

We have limited our study to models with 7 = 0, 1/2, 1, 3/2, and 2, which represent a wide 
range for the density cusp at small radii. The 7 = model is the only one of this family that has 
a (non-isothermal) core. The model that most closely resembles the R^^'^ law in surface density is 
the 7 = 3/2 model (D93). In spite of a previous study by Merritt & Aguilar (1985) of JafFe (7 = 2) 
model, we have included it here to obtain a more accurate estimation for its stability threshold. 
In fact, as we show below, we found a slightly larger value for the critical anisotropy radius than 
theirs. 

For all these models, we have generated a set of initial conditions using the Osipkov-Merritt 
DF|^. Each model was truncated at a radius enclosing 99.9% of the total mass. All simulations 
employ N = 50, 000 equal mass particles. Table |^ shows the parameters associated with each of 
these models; where rh is the half-mass radius, rg is the lower value for the anisotropy radius such 
that the DF is strictly positive, and Th is the dynamical time evaluated at rh- 



2.2. N-body code 

The stability of these models has been studied using an N-body code based on the 'self- 
consistent field' method developed by H092 (see also Clutton-Brock 1973). In this approach, the 
density and the gravitational potential are expanded in a biorthogonal set of basis functions as 

Pi^) = ^^nlmPnlm{Tc) = ^ AnlmPnl{r)Yimi9 , if) , 
nlm nlm 

(3) 

$(r) = ^Anlm^nlmir) = ^ ^„,Z,ri^ni ('^)^m V') ■ 
nlm nlm 

II092 obtained their basis functions {pnimi ^nim} using the model for spherical galaxies pro- 
posed by II90 as the zeroth-order term. Higher order terms are found by construction. For a 
set of N pointlike particles the expansion coefficients Anim are given in terms of the positions of 
the particles. Then, the accelerations are computed by a simple analytical differentiation of the 
potential expansion (see II092 for a more detailed description). 

One of us (AM) has written a code that implements this formalism. In this code the particle 
positions and velocities are updated using a second order integrator with a fixed timestep At, given 



There is a factor 2 missing in the second term of the right-hand side of equation (A7) of D93. 
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by 

Xi+i = Xi + AtVi + ^At^ai, (4) 

Vj+i = Vj + ^ At (aj + a^+i), (5) 

where the subscript identifies the iteration. This integrator is equivalent to the standard time 
centered leapfrog, as can be verified by direct substitution (Allen & Tildesley 1992; Hut, Makino, 

6 McMillan 1995). We have studied its numerical properties and found them closely similar to 
those of the standard leapfrog. 

The H092 basis set can be used to represent a number of potential-density pairs for spherical 
galaxies, including the 7-models. For example, the 7 = 1 model is just the zeroth order term and the 

7 = model can be reproduced exactly as a linear combination of two terms of this basis (H092). 
However, for practical purposes, the potential (density) of other 7-models are approximated using 
truncated expansions of the kind defined in equation dH) limiting the number of radial functions to 

In our simulations, nmax has been chosen according to the studied model. Figure || shows 
the values of the expansion coefficients A^qo (see eq. [2.35] of H092) for some of the 7-models. 
The coefficients are directly related with 7; smaller values of Anoo are obtained for lower 7. In all 
the cases, for large values of n the coefficients Anoo levels off. Then, adding more terms to the 
potential expansion does not decrease significantly the error, but the computational time increases 
considerably. However, in most cases, good accuracy in the evaluation of the potential can be 
achieved with rimax ~ 6. This point is illustrated in Figure ^, where we compare the radial 
acceleration computed using the H092 basis functions for the potential expansion (^), truncated 
at different values of nmaxj to its exact value for the 7 = 1/2 and 7 = 3/2 models. In both cases, at 
large radii, the acceleration obtained using the truncated potential expansion agrees with the one 
obtained using the full potential. However, at small radii, the error grows with 7. For 7 = 1/2, 
the radial acceleration can be computed to better than 0.5% accuracy at small radii, while for 7 = 
3/2 the same quantity can be evaluated only up to 5% accuracy. The adopted values for rimax (see 
Table |) were chosen such that the error in the acceleration, at small radii, was lower that 5%. 

Since we are interested in the ROI, all our simulations only include series expansions up to 
^max = 2. The total elapsed time was chosen as 50 half-mass dynamical times. The timestep for the 
different models appears in Table |^. They were taken such that the energy was conserved better 
than 0.5% during the total elapsed time. In some special cases, we checked that our results would 
not depend on the number of terms used in the potential expansion and the timestep (details are 
given below). 
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3. Results 



To study the ROI in the models introduced in § 2.1, we have done a set of simulations with 
different values for the anisotropy radius ra- For each run, we have used an iterative algorithm to 
estimate the axial ratios of the particle distribution. In this scheme, initial values for the modified 
inertia tensor 

^^.=E5^' a' = x'+y'/<lf + zVqi (6) 

are calculated for all particles inside a spherical shell (^i = = 1) with a given radius r^. New 
axis ratios qi, q2 and the orientation of the fitting ellipsoid are then estimated from the eigenvalues 
of lij and used to obtain an improved approximation to the modified inertia tensor. This process 
is repeated until the axis ratios converge to a value within a pre-established tolerance, Aq = 0.001. 

We have tested the precision of this scheme computing the axis ratio of a set of random N 
particle positions generated from the following triaxial generalization of the Hernquist model (e.g., 
Merritt & Fridman 1996) 

^("^^ = 2^^(1 + ^)3' 
with ^ ^ 

= + ^ + ^, 0<c<6<l. (8) 

Figure ^ shows the average computed minor axis ratio as a function of the number of particles 
within the radius r^. The error bar quoted correspond to the standard deviation for 10 different 
measures. We found that for ^ 1000, the computed axis ratio of the particle distribution came 
out with an error smaller than 1%. Convergence problems appear only when the number of particles 
inside the measuring radius is small, or when itself is small. To avoid this situation we have 
adopted as a measuring radius the one that encloses the 70% of the total mass of the model. 
With these values for r^, we have obtained the axis ratio used as a test for stability. 

The evolution of the ratio between the minor to major axis for the Hernquist model with 
different values of the anisotropy radius appears on Figure These values were obtained by 
fitting an ellipsoid at radius = 5. The stability boundary appears to be close to = 1.0 
and we assume that this model is stable. Models with ^ 0.8 are strongly unstable and their 
final configuration is nearly prolate, with the final axis ratio determined by their initial anisotropy; 
the larger the value of the initial anisotropy, the smaller the final axis ratios. An example of the 
final state of a strongly unstable Hernquist model appears on Figure |5|, where we have plotted 
the particle positions, projected along the minor semiaxis (z-axis), at the start and the end of a 
simulation for a model with anisotropy radius = 0.3. The particle positions are rotated in such 
a way that the semiaxis become aligned with the cartesian axis. At the end of the run, it is clearly 
visible a (triaxial) bar at the center of the system with c/a = 0.33 and b/a = 0.40. 

Figure ^ shows the evolution of the minor to major axis ratio for the Jaffe (7 = 2) model. 
The axis ratios were evaluated at radius = 2.3. In their study of this model, Merritt &; Aguilar 
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(1985) found a critical value for between 0.2 and 0.3. However, the total elapsed time in their 
simulations was only 20 dynamical times. This ambiguity is settled in Figure ^, where we can see 
that the model with = 0.3 appears to be stable for elapsed times t ^ 20Th, however, for t ~ SOT^ 
the system show signs of an incipient bar. Accordingly, we set the stability threshold at 0.4. 

The evolution of the axis ratio for other 7-models is very similar. Our results are summarized 
in Figure ^. Table |l| includes The critical anisotropy radius and the respective critical global 
anisotropy 2Tr/Tt. The stability boundary between stable and unstable models is given roughly by 
the linear fit = 1.54 — O.547. This boundary came out very close with the conservative estimate 
Va = 1.6 — 0.677, raised by Carollo, de Zeeuw, & van der Marel (1995) based on the supposition 
that the stability boundary for models with 7 < 2 is given by 2Tr/Tt = 2.5, the critical global 
anisotropy quoted by Merritt & Aguilar (1985) for the 7 = 2 model. 

Our results give a critical global anisotropy parameter 2Tr/T± = 2.31 =b 0.27 for the 7-models; 
however, the values of the anisotropy parameter are larger for lower values of 7. This limit is larger 
than the values 1.75 and 1.58 reported by Fridman &: Polyachenko (1984) and Bertin et al. (1994), 
respectively. In our opinion, these results strengthen the belief that this criterion does not yield a 
reliable prediction for the stability threshold value because it is strongly model dependent. 

We did additional runs to verify that our results were indeed independent of the timestep and 
the number of terms used in the potential expansion. Figure ^ shows the evolution of the minor to 
major axis ratio for the Hernquist (7 = 1) model with anisotropy radius = 0.3. The run for this 
model was repeated using a shorter timestep. At = 0.01, with essentially the same behavior for 
the axis ratio. A run with rimax = 6 and /max = 4 converged to a slightly smaller final axis ratio. 
Similar results were obtained for other examples. 



3.1. Radial stability 

The sufficient condition df/dE < (Doremus & Feix 1973; Dejonghe & Merritt 1988) can 
be used to demonstrate the stability against radial perturbations of the Osipkov- Merritt 7-models. 
The lower value for the anisotropy radius which satisfies the Doremus-Feix's criterion, r^p, is given 
in Table |l[ We have studied the stability to radial perturbation modes of some of those models 



that fail to satisfy this criterion, using the same N-body code described in § 2.2, but retaining only 
the I = terms in the potential expansion. No signs of radial instability were observed in these 
runs. 

An example of this result is displayed in Figure P, which is a plot of the evolution of the 
radii containing 10%, 20%, . . . , 70% of the total mass and the minor to major axis ratio, for the 
7 = model with anisotropy radius = 0.5. In both cases, the mass distribution only shows 
the fluctuations associated with the finite number of particles. There are no discernible changes in 
the radial distribution of matter and the shape of the particle distribution. Similar results were 
observed in the other examples, which do not satisfy the Doremus-Feix's criterion. We conclude 
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that all these models are stable to radial modes. 



4. Discussion and summary 

Recently, Hjorth (1994) has proposed an analytical criterion for the onset of instability for 
Osipkov-Merritt models. We have used our simulations to test this criterion. For the models 
discussed here, we have found that the numerical instability thresholds are higher than the predicted 
value given by Hjorth's criterion, df /dQ < (circles in Figure |^. Other independent numerical 
results point in the same direction. Dejonghe & Merritt (1988) using a multipolar expansion 
N-body code have found that the Osipkov-Merritt-Plummer model is unstable for an anisotropy 
radius ^ 1.1, while Hjorth's criterion predicts ^ 0.9. For this model, we also have tried a 
different approach doing a series of simulations using an N-body code based on the Glutton-Brock 
(1973) basis set (see also H092). Our results (Meza &: Zamorano 1996) indicate that the stability 
boundary lies near — 1.1? the same value previously obtained by Dejonghe Merritt (1988). 

According to Hjorth (1994), this discrepancy may be originated by numerical errors in the 
representation of the central potential due to the small number of particles used in the simula- 
tions. He argues that these errors could induce artificial inflection points in the DF of the N-body 
realization. In this case, there should be a correlation between the number of particles and the 
stability threshold determined using these conditions: fewer particles should start an instability for 
higher values of the anisotropy radius r^. To decide about this eventual behavior we have done 
some additional runs with N = 10^ and 2 x 10^ particles for anisotropy radius closer to the critical 



value quoted for the 7-models studied in the previous section. Figure |10| shows the results for the 
evolution of the minor to major axis ratio for the 7 = 3/2 model with anisotropy radius = 0.8. 
We observe that the curves show basically the same behavior and we conclude that the number 
of particles does not affect our estimation for the stability threshold. The same result repeats for 
the other 7-models. After these considerations, we feel that it would be interesting to review the 
hypothesis used in the Hjorth's criterion. 

Our results can be summarized as follows: 



1. Using an N-body code based on the 'self-consistent field' method developed by H092, we have 
tested the stability of a set of 7-models with DFs of the Osipkov-Merritt type. We found an 
approximated stability boundary for the ROI on the (7, ra)-plane. This boundary is given by 
a global anisotropy parameter 2Tr/Tt = 2.31. 

2. The criterion df /dQ < 0, recently suggested by Hjorth, gives a lower stability threshold for 
the 7-models. 

3. We have studied the stability to radial modes of those 7-models which fail to satisfy the 
Doremus-Feix's stability criterion. No signs of instability were observed and we conclude that 
these models are stable to radial (shell forming) modes. 
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Table 1. Parameters of the models 
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Fig. 1.— 



Expansion coefficients for some 7-models using tfie HO 92 basis set. 




Fig. 2. — Relative error in the radial acceleration for the 7 = 1/2 (a) and 7 = 3/2 (b) model as 
a function of radius r. The figures compare the acceleration, as computed using the H092 basis 
functions for different values of nmaxi to their exact values. 
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Fig. 3. — Computed minor axis ratio as a function of the number of particles N for a set of 10 
random generations of the "triaxial" Hernquist profile (eq. [^) with c = 0.3 (down triangles), 
c = 0.5 (up triangles), c = 0.7 (circles), and c = 0.9 (squares). The error bar quoted here 
corresponds to the standard deviation of the measures. 
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Fig. 4. — Evolution of the minor to major axis ratio for the Hernquist (7 = 1) model for different 
values of the anisotropy radius Tq. The axis ratios were obtained by fitting an ellipsoid at radius 
= 5, the radius which encloses ~ 70% of the total mass. The time is normalized to the half-mass 
dynamical time T^. Solid line: = 1.1; long-dashed line: = 1.0; short-dashed line: Va = 0.9; 
dotted-dashed line: = 0.7; and dotted line: = 0.3. 
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Fig. 5. — Particle positions as viewed along the minor semiaxis at the beginning (a) and at the 
end (b) of a simulation for the Hernquist model with Va = 0.3. The particle positions are rotated 
such that the major (intermediate) semiaxis remains aligned with the x-axis (y-axis). This region 
contains approximately 70% of the total mass of the system. 
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Fig. 6. — Evolution of the minor to major axis ratio for the Jaffe (7 = 2) model. These values were 
obtained by fitting an ellipsoid at radius rm = 2.3, the radius which encloses ~ 70% of the total 
mass. The time is normalized to the half-mass dynamical time T^. Solid line: = 0.4; dashed 
line: Va = 0.3; and dotted line: Va = 0.2. 
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Fig. 7. — The (7, ra)-plane for Osipkov-Merritt 7-models. Models in the shaded region have 
f{Q) < for some values of Q. The squares plotted are our numerical results and the stability 
threshold predicted by the Hjorth's criterion df/dQ < is marked with circles. The solid curve 
drawn is a linear fit obtained with our numerical results and roughly represents the lower boundary 
for models which are not affected by the radial-orbit instability. 
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Fig. 8. — Evolution of the minor to major axis ratio for a Hernquist (7 = 1) model with = 0.3. 
The same quantity is given for simulations with a shorter timestep At = 0.01, and a potential 
expansion with rimax = 6 and Zmax = 4. 
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Fig. 9. — Evolution of the radial distribution of matter (a) and the axis ratios (b) for the 7 = 
model with anisotropy radius = 0.5. These results were obtained retaining only the I = terms 
in the potential expansion. 
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Fig. 10. — Evolution of the minor to major axis ratio for the 7 = 3/2 model using N = 5x 10*^, 10^, 
and 2 x 10^ particles and anisotropy radius Va = 0.8. These values were obtained fitting an ellipsoid 
at radius rm = 3.7. 



